<!DOCTYPE HTML PUBLIC "-//W3C//DTD HTML 4.01 Transitional//EN"
                "http://www.w3.org/TR/REC-html40/loose.dtd">
<html>
<head>
  <title>Description of v_imagehomog</title>
  <meta name="keywords" content="v_imagehomog">
  <meta name="description" content="V_IMAGEHOMOG Apply a homography transformation to an image with bilinear interpolation">
  <meta http-equiv="Content-Type" content="text/html; charset=iso-8859-1">
  <meta name="generator" content="m2html &copy; 2003 Guillaume Flandin">
  <meta name="robots" content="index, follow">
  <link type="text/css" rel="stylesheet" href="../m2html.css">
</head>
<body>
<a name="_top"></a>
<div><a href="../index.html">Home</a> &gt;  <a href="index.html">v_mfiles</a> &gt; v_imagehomog.m</div>

<!--<table width="100%"><tr><td align="left"><a href="../index.html"><img alt="<" border="0" src="../left.png">&nbsp;Master index</a></td>
<td align="right"><a href="index.html">Index for v_mfiles&nbsp;<img alt=">" border="0" src="../right.png"></a></td></tr></table>-->

<h1>v_imagehomog
</h1>

<h2><a name="_name"></a>PURPOSE <a href="#_top"><img alt="^" border="0" src="../up.png"></a></h2>
<div class="box"><strong>V_IMAGEHOMOG Apply a homography transformation to an image with bilinear interpolation</strong></div>

<h2><a name="_synopsis"></a>SYNOPSIS <a href="#_top"><img alt="^" border="0" src="../up.png"></a></h2>
<div class="box"><strong>function [ih,xa,ya]=v_imagehomog(im,h,m,clip) </strong></div>

<h2><a name="_description"></a>DESCRIPTION <a href="#_top"><img alt="^" border="0" src="../up.png"></a></h2>
<div class="fragment"><pre class="comment">V_IMAGEHOMOG Apply a homography transformation to an image with bilinear interpolation
Inputs: im(ny,nx,nc)  input image (uint8)
        h(3,3)        homography
        m             mode string
                         i  show image [default if no output arguments]
                         m  h uses matlab coordinates: (1,1) is top left [default]
                         c  h uses central coordinates: (0,0) is the centre = {(1+nx)/2,(1+ny)/2} in 'm'
                         k  clip to original image dimensions (equivalent to clip=1)
                         x  add blank rows and cols to extend image to the specified clipping rectangle
                         t  trim output image by re moving blank rows and columns
        clip(4)       output image bounding box [xmin xmax ymin ymax] (same coordinate system as h)
                         or [n] or [nx ny] to clip to a multiple of the original image [default clip=2]
 Outputs:
        ih(my,mx,nc)  output image (uint8)
        xa(mx)        x axis
        ya(my)        y axis</pre></div>

<!-- crossreference -->
<h2><a name="_cross"></a>CROSS-REFERENCE INFORMATION <a href="#_top"><img alt="^" border="0" src="../up.png"></a></h2>
This function calls:
<ul style="list-style-image:url(../matlabicon.gif)">
</ul>
This function is called by:
<ul style="list-style-image:url(../matlabicon.gif)">
<li><a href="v_rectifyhomog.html" class="code" title="function [imr,xa,ya]=v_rectifyhomog(ims,roc,k0,mode)">v_rectifyhomog</a>	V_RECTIFYHOMOG Apply rectifying homographies to an image set</li></ul>
<!-- crossreference -->


<h2><a name="_source"></a>SOURCE CODE <a href="#_top"><img alt="^" border="0" src="../up.png"></a></h2>
<div class="fragment"><pre>0001 <a name="_sub0" href="#_subfunctions" class="code">function [ih,xa,ya]=v_imagehomog(im,h,m,clip)</a>
0002 <span class="comment">%V_IMAGEHOMOG Apply a homography transformation to an image with bilinear interpolation</span>
0003 <span class="comment">%Inputs: im(ny,nx,nc)  input image (uint8)</span>
0004 <span class="comment">%        h(3,3)        homography</span>
0005 <span class="comment">%        m             mode string</span>
0006 <span class="comment">%                         i  show image [default if no output arguments]</span>
0007 <span class="comment">%                         m  h uses matlab coordinates: (1,1) is top left [default]</span>
0008 <span class="comment">%                         c  h uses central coordinates: (0,0) is the centre = {(1+nx)/2,(1+ny)/2} in 'm'</span>
0009 <span class="comment">%                         k  clip to original image dimensions (equivalent to clip=1)</span>
0010 <span class="comment">%                         x  add blank rows and cols to extend image to the specified clipping rectangle</span>
0011 <span class="comment">%                         t  trim output image by re moving blank rows and columns</span>
0012 <span class="comment">%        clip(4)       output image bounding box [xmin xmax ymin ymax] (same coordinate system as h)</span>
0013 <span class="comment">%                         or [n] or [nx ny] to clip to a multiple of the original image [default clip=2]</span>
0014 <span class="comment">% Outputs:</span>
0015 <span class="comment">%        ih(my,mx,nc)  output image (uint8)</span>
0016 <span class="comment">%        xa(mx)        x axis</span>
0017 <span class="comment">%        ya(my)        y axis</span>
0018 
0019 <span class="comment">% Bugs/Suggestions:</span>
0020 <span class="comment">% (1) cope with non-uint8</span>
0021 <span class="comment">% (2) cope with (a) multiple inputs, (b) multiple transformations</span>
0022 <span class="comment">% (3) output a boundary mask as an alpha channel</span>
0023 <span class="comment">% (4) do anti-aliasing along the boundary</span>
0024 <span class="comment">% (5) check that origin shift is correct for central coordinates</span>
0025 
0026 <span class="comment">%      Copyright (C) Mike Brookes 2010-2012</span>
0027 <span class="comment">%      Version: $Id: v_imagehomog.m 10865 2018-09-21 17:22:45Z dmb $</span>
0028 <span class="comment">%</span>
0029 <span class="comment">%   VOICEBOX is a MATLAB toolbox for speech processing.</span>
0030 <span class="comment">%   Home page: http://www.ee.ic.ac.uk/hp/staff/dmb/voicebox/voicebox.html</span>
0031 <span class="comment">%</span>
0032 <span class="comment">%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%</span>
0033 <span class="comment">%   This program is free software; you can redistribute it and/or modify</span>
0034 <span class="comment">%   it under the terms of the GNU General Public License as published by</span>
0035 <span class="comment">%   the Free Software Foundation; either version 2 of the License, or</span>
0036 <span class="comment">%   (at your option) any later version.</span>
0037 <span class="comment">%</span>
0038 <span class="comment">%   This program is distributed in the hope that it will be useful,</span>
0039 <span class="comment">%   but WITHOUT ANY WARRANTY; without even the implied warranty of</span>
0040 <span class="comment">%   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the</span>
0041 <span class="comment">%   GNU General Public License for more details.</span>
0042 <span class="comment">%</span>
0043 <span class="comment">%   You can obtain a copy of the GNU General Public License from</span>
0044 <span class="comment">%   http://www.gnu.org/copyleft/gpl.html or by writing to</span>
0045 <span class="comment">%   Free Software Foundation, Inc.,675 Mass Ave, Cambridge, MA 02139, USA.</span>
0046 <span class="comment">%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%</span>
0047 
0048 maxby=1e7;  <span class="comment">% maximum memory to use</span>
0049 [ny,nx,nc]=size(im);
0050 <span class="keyword">if</span> nargin&lt;4 || ~numel(clip)
0051     clip=2;
0052 <span class="keyword">end</span>
0053 <span class="keyword">if</span> nargin&lt;3 || ~numel(m)
0054     m=<span class="string">''</span>;
0055 <span class="keyword">end</span>
0056 <span class="keyword">if</span> nargin&lt;2 || ~numel(h)
0057     h=eye(3);
0058 <span class="keyword">end</span>
0059 
0060 imr=reshape(im,nx*ny,nc); <span class="comment">% vectorize the input image</span>
0061 t=eye(3);
0062 <span class="keyword">if</span> any(m==<span class="string">'c'</span>)   <span class="comment">% convert homography and clipping box to matlab coordinates</span>
0063     t(7:8)=0.5+[nx ny]/2;  <span class="comment">% shift origin to image centre</span>
0064     h=t*h/t;  <span class="comment">% change homography so input and output use MATLAB coordinates</span>
0065     <span class="keyword">if</span> numel(clip)==4
0066         clip=clip+t([7 7 8 8]);  <span class="comment">% make clipping use MATLAB coordinates as well</span>
0067     <span class="keyword">end</span>
0068 <span class="keyword">end</span>
0069 <span class="comment">% determine clipping rectangle for output image</span>
0070 <span class="keyword">if</span> any(m==<span class="string">'k'</span>)
0071     clip=[1 1];
0072 <span class="keyword">elseif</span> numel(clip)==1
0073     clip=clip(ones(1,2));
0074 <span class="keyword">end</span>
0075 <span class="keyword">if</span> numel(clip)==2
0076     clip=[[1 nx]*clip(1)-(clip(1)-1)*(1+nx)/2 [1 ny]*clip(2)-(clip(2)-1)*(1+ny)/2];
0077 <span class="keyword">end</span>
0078 clip=clip(:)';
0079 clip(1:2:3)=floor(clip(1:2:3));
0080 clip(2:2:4)=ceil(clip(2:2:4));
0081 <span class="comment">% now determine the image of the source image</span>
0082 bi=[1 1 nx nx; 1 ny ny 1; 1 1 1 1];
0083 box=h*bi;
0084 b3=box(3,:);
0085 <span class="keyword">if</span> any(b3&lt;=0)
0086     ib=find(b3&gt;0);
0087     nb=numel(ib);
0088     bb=ones(3,nb+2);
0089     <span class="keyword">if</span> ~nb
0090         error(<span class="string">'image invisible'</span>);
0091     <span class="keyword">end</span>
0092     bb(:,1:nb)=box(:,ib);
0093     px=[4 1:3];
0094     ib=find(b3.*b3(px)&lt;=0); <span class="comment">% find points at infinity</span>
0095     ip=px(ib);
0096     af=b3(ip);
0097     bf=b3(ib);
0098     pof=[-1 0 1 0; 0 1 0 -1; 0 0 0 0]; <span class="comment">% offsets to avoid exact infinity</span>
0099     w3=ones(3,1);
0100     bb(:,nb+(1:2))=h*((bf(w3,:).*bi(:,ip)-af(w3,:).*bi(:,ib))./repmat(bf-af,3,1)+(pof(:,ib).*repmat(2*(b3(ib)&gt;0)-1,3,1)));
0101     box=bb;
0102 <span class="keyword">end</span>
0103 box=box(1:2,:)./box([3 3],:);  <span class="comment">% coordinates of mapped original image</span>
0104 box=[min(box(1,:)) max(box(1,:)) min(box(2,:)) max(box(2,:))];
0105 
0106 box(1:2:3)=floor(max(clip(1:2:3),box(1:2:3)));  <span class="comment">% no point in calculating non-existant points</span>
0107 box(2:2:4)=ceil(min(clip(2:2:4),box(2:2:4)));
0108 g=inv(h);
0109 mx=box(2)-box(1)+1; <span class="comment">% number of columns in destination</span>
0110 my=box(4)-box(3)+1; <span class="comment">% number of rows in destination</span>
0111 
0112 ih=zeros(my*mx,nc,<span class="string">'uint8'</span>);
0113 ncol=max(1,min(mx,floor(maxby/(my*nc*8)))); <span class="comment">% number of columns to do in a chunk</span>
0114 nloop=ceil(mx/ncol);
0115 cmax=1+rem(mx-1,ncol); <span class="comment">% final column of first iteration</span>
0116 jxinc=ncol*my;      <span class="comment">% increment target indices each loop</span>
0117 ginc=g(:,1)*ncol; <span class="comment">% increment transformed targets each loop</span>
0118 wj=ones(1,jxinc);   <span class="comment">% repeat index for ginc</span>
0119 jx=1+cmax*my-jxinc:cmax*my;  <span class="comment">% initial target indices (some might be negative)</span>
0120 gk=g*[reshape(repmat(cmax-ncol+box(1):cmax+box(1)-1,my,1),1,jxinc); repmat(box(3):box(4),1,ncol); ones(1,jxinc)];
0121 gn=gk(1:2,:)./gk([3 3],:);   <span class="comment">% normalize source coordinates</span>
0122 mn=[zeros(1,jxinc-cmax*my) ones(1,cmax*my)]; <span class="comment">% mask for initial iteration</span>
0123 mn=mn &amp; (gn(1,:)&gt;-0.5 &amp; gn(2,:)&gt;-0.5 &amp; gn(1,:)&lt;nx+0.5 &amp; gn(2,:)&lt;ny+0.5); <span class="comment">% mask valid pixels</span>
0124 w3=ones(nc,1);
0125 <span class="keyword">for</span> i=1:nloop
0126     fn=max(floor(gn(:,mn)),1);
0127     fn1=min(max(fn(1,:)',1),nx-1);
0128     fn2=min(max(fn(2,:)',1),ny-1);
0129     dn=gn(:,mn)-[fn1 fn2]';
0130     dn1=min(max(dn(1,:)',0),1);
0131     dn2=min(max(dn(2,:)',0),1);
0132     dn1c=1-dn1;
0133     dn2c=1-dn2;
0134     ih(jx(mn),:)=uint8(dn1c(:,w3).*(dn2c(:,w3).*single(imr(fn2+ny*(fn1-1),:)) <span class="keyword">...</span>
0135         +dn2(:,w3).*single(imr(fn2+1+ny*(fn1-1),:))) <span class="keyword">...</span>
0136         +dn1(:,w3).*(dn2c(:,w3).*single(imr(fn2+ny*(fn1),:)) <span class="keyword">...</span>
0137         +dn2(:,w3).*single(imr(fn2+1+ny*(fn1),:))));
0138     jx=jx+jxinc;  <span class="comment">% target indices</span>
0139     gk=gk+ginc(:,wj);
0140     gn=gk(1:2,:)./gk([3 3],:);   <span class="comment">% normalize source coordinates</span>
0141     mn=gn(1,:)&gt;-0.5 &amp; gn(2,:)&gt;-0.5 &amp; gn(1,:)&lt;nx+0.5 &amp; gn(2,:)&lt;ny+0.5; <span class="comment">% mask valid pixels</span>
0142 <span class="keyword">end</span>
0143 ih=reshape(ih,[my,mx,nc]);
0144 <span class="keyword">if</span> any(m==<span class="string">'x'</span>)          <span class="comment">% extend blank area to specified clipping rectangle</span>
0145     ih=[zeros(box(3)-clip(3),clip(2)-clip(1)+1,nc,<span class="string">'uint8'</span>); <span class="keyword">...</span>
0146         zeros(my,box(1)-clip(1),nc,<span class="string">'uint8'</span>) ih zeros(my,clip(2)-box(2),nc,<span class="string">'uint8'</span>);
0147         zeros(clip(4)-box(4),clip(2)-clip(1)+1,nc,<span class="string">'uint8'</span>)];
0148     xa=(clip(1):clip(2))-t(7);
0149     ya=(clip(3):clip(4))-t(8);
0150 <span class="keyword">else</span>
0151     xa=(box(1):box(2))-t(7);
0152     ya=(box(3):box(4))-t(8);
0153 <span class="keyword">end</span>
0154 <span class="keyword">if</span> any(m==<span class="string">'t'</span>) 
0155      ihs=sum(ih,3);
0156      azx=all(~ihs,1);
0157      ix1=find(~azx,1,<span class="string">'first'</span>);
0158      <span class="keyword">if</span> ~numel(ix1)
0159          ih=[];
0160          xa=[];
0161          ya=[];
0162      <span class="keyword">else</span>
0163          ix2=find(~azx,1,<span class="string">'last'</span>);
0164          azx=all(~ihs,2);
0165          iy1=find(~azx,1,<span class="string">'first'</span>);
0166          iy2=find(~azx,1,<span class="string">'last'</span>);
0167          ih=ih(iy1:iy2,ix1:ix2,:);
0168          xa=xa(ix1:ix2);
0169          ya=ya(iy1:iy2);
0170      <span class="keyword">end</span>
0171 <span class="keyword">end</span>
0172 <span class="keyword">if</span> ~nargout || any(m==<span class="string">'i'</span>)
0173     imagesc(xa,ya,ih);
0174     axis image
0175 <span class="keyword">end</span>
0176</pre></div>
<hr><address>Generated by <strong><a href="http://www.artefact.tk/software/matlab/m2html/">m2html</a></strong> &copy; 2003</address>
</body>
</html>